library(QMSS)
d = read.csv(file.choose())
library(plyr)
library(leaps)

library(stargazer)

vars <- c("idnum","panelwave",
          "happy",
          "age",
          "health",
          "marital",
          "income",
          "realinc",
          "sex",
          "race",
          "degree",
          "richwork",
          "finalter",
          "getahead",
          "parsol")
pd.sub <- d[, vars]

#recode happy
pd.sub$rhappy[pd.sub$happy==1] <- 3
pd.sub$rhappy[pd.sub$happy==2] <- 2
pd.sub$rhappy[pd.sub$happy==3] <- 1
#recode race
pd.sub$white =ifelse(pd.sub$race==1,1,0)
#recode health
pd.sub$rhealth[pd.sub$health==1]=4
pd.sub$rhealth[pd.sub$health==2]=3
pd.sub$rhealth[pd.sub$health==3]=2
pd.sub$rhealth[pd.sub$health==4]=1
#recode marital 
pd.sub$brokenmar[pd.sub$marital==2]=1
pd.sub$brokenmar[pd.sub$marital==3]=1
pd.sub$brokenmar[pd.sub$marital==4]=1
pd.sub$brokenmar[pd.sub$marital==1]=0
pd.sub$brokenmar[pd.sub$marital==5]=0
#recode richwork
pd.sub$workformoney[pd.sub$richwork==1]=0
pd.sub$workformoney[pd.sub$richwork==2]=1

#recode getahead
pd.sub$hardwork=ifelse(pd.sub$getahead==1,1,0)
#recode finalter
pd.sub$finalterbetter=ifelse(pd.sub$finalter==1,1,0)

#recode parsol
pd.sub$parsolbetter=ifelse(pd.sub$parsol==1 | pd.sub$parsol==2,1,0)
#recode wave
pd.sub$year[pd.sub$panelwave==1] = 2006 
pd.sub$year[pd.sub$panelwave==2] = 2008
pd.sub$year[pd.sub$panelwave==3] = 2010


library(psych)
describe(pd.sub, skew = FALSE )
Tab(pd.sub$parsolbetter)
Tab(pd.sub$workformoney)
Tab(pd.sub$hardwork)

library(ggplot2)

ggplot(d, aes(x = age))+
  geom_line(stat = "density", colour = 'red')+theme(axis.title.x = element_text(size=40,colour = "black"),axis.title.y = element_text(size=40,colour = "black"),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank())+theme(axis.text=element_text(size=20))

ggplot(d, aes(x = realinc))+
  geom_line(stat = "density", colour = 'red')+theme(axis.title.x = element_text(size=40,colour = "black"),axis.title.y = element_text(size=40,colour = "black"),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank())+theme(axis.text=element_text(size=20))+
  labs(x='Real Income')


ggplot(d, aes(x = factor(income)))+
  geom_histogram(fill = 'white', colour='black')+theme(axis.title.x = element_text(size=40,colour = "black"),axis.title.y = element_text(size=40,colour = "black"),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank())+theme(axis.text=element_text(size=20))+
  labs(x = "Income Level")


ggplot(pd.sub, aes(x = rhappy)) + geom_histogram(fill = 'white', colour='black', position = 'dodge')+
  labs(x = "Happiness")+
  theme(axis.title.x = element_text(size=50,colour = "black"),axis.title.y = element_text(size=50,colour = "black"),panel.grid.major = element_blank())+
  theme(axis.text=element_text(size=20))+
  theme(strip.text = element_text(face = "bold", size=rel(3.5)))+
  facet_grid(year~.)



pd.sub$f[pd.sub$finalterbetter==1]='Better'
pd.sub$f[pd.sub$finalterbetter==0]='Not Better'
pd.sub$happiness = as.factor(pd.sub$rhappy)
ggplot(pd.sub, aes(x = happiness)) + geom_histogram(fill = 'white', colour='black', position = 'dodge')+
  labs(x = "Happiness", titile = 'Financial Situation')+
  theme(axis.title.x = element_text(size=30,colour = "black"),axis.title.y = element_text(size=30,colour = "black"),panel.grid.major = element_blank())+
  theme(axis.text=element_text(size=20))+
  theme(strip.text = element_text(face = "bold", size=rel(3.5)))+
  facet_grid(.~f)

ggplot(pd.sub, aes(x = rhappy)) + geom_histogram(fill = 'white', colour='black', position = 'dodge')+
  labs(x = "Happiness")+
  theme(axis.title.x = element_text(size=50,colour = "black"),axis.title.y = element_text(size=50,colour = "black"),panel.grid.major = element_blank())+
  theme(axis.text=element_text(size=20))+
  theme(strip.text = element_text(face = "bold", size=rel(3.5)))+
  facet_grid(.~hardwork)




lm0<-lm(rhappy~workformoney+parsolbetter+hardwork+rhealth+finalterbetter+log(realinc)+age+brokenmar+sex+white+degree,data = pd.sub)
summary(lm0)
stargazer(lm0, type = "text")

plm1<-plm(rhappy~workformoney+parsolbetter+hardwork+rhealth+finalterbetter+log(realinc)+age+brokenmar+sex+white+degree,index = c("idnum", "panelwave"),  model = "pooling",data = pd.sub)
summary(plm1)
clusterSE(fit = plm1, cluster.var = "idnum", data=pd.sub)




pd.sub$d.age<-firstD(pd.sub$age)
pd.sub$d.sex<-firstD(pd.sub$sex)
pd.sub$d.white<-firstD(pd.sub$white)

pd.sub$d.rhappy<-firstD(pd.sub$rhappy)
pd.sub$d.workformoney<-firstD(pd.sub$workformoney)
pd.sub$d.parsolbetter<-firstD(pd.sub$parsolbetter)
pd.sub$d.hardwork<-firstD(pd.sub$hardwork)
pd.sub$d.rhealth<-firstD(pd.sub$rhealth)
pd.sub$d.income<-firstD(pd.sub$realinc)
pd.sub$d.finalterbetter<-firstD(pd.sub$finalterbetter)
pd.sub$d.brokenmar<-firstD(pd.sub$brokenmar)
pd.sub$d.degree<-firstD(pd.sub$degree)






#stepwise
stp <- regsubsets(d.rhappy~d.workformoney+d.parsolbetter+d.hardwork+d.rhealth+d.finalterbetter+d.income+d.age+d.brokenmar+d.sex+d.white+d.degree,data = pd.sub,nbest = 1, nvmax = 10)
plot(stp, scale="bic")
plot(stp, scale="adjr2")

ran <- plm(formula = rhappy~workformoney+parsolbetter+hardwork+rhealth+log(realinc)+brokenmar+white,
           index = c("idnum", "panelwave"),  model = "random",
           data = pd.sub)
summary(ran)
fix <- plm(formula = rhappy~workformoney+parsolbetter+hardwork+rhealth+log(realinc)+brokenmar+white,
           index = c("idnum", "panelwave"),  model = "within",
           data = pd.sub)
summary(fix)
Hausman<-phtest(fix, ran)

Hausman
plm3<-plm(rhappy~workformoney+parsolbetter+hardwork+rhealth+log(realinc)+brokenmar+white,
          index = c("idnum", "panelwave"),  model = "fd",
          data = pd.sub)
summary(plm3)

stargazer(fix, ran, plm3, type = 'text')

